﻿#if 1
#include <math.h>
#include "novaFractal.h"
#include "../common/glBitmap.h"
#include "../common/common.h"

/*==========================================================================*/

/*==========================================================================*/

struct nvComplex
{
	float r;
	float i;

	nvComplex(float a, float b) : r(a), i(b) 
	{ }

	__device__ float squareLength() 
	{
		return r * r + i * i;
	}

	__device__ nvComplex operator*(const float a)
	{
		return nvComplex(r*a, i*a);
	}
	__device__ nvComplex operator*(const nvComplex& a)
	{
		return nvComplex(r*a.r - i*a.i, i*a.r + r*a.i);
	}

	__device__ nvComplex operator+(const nvComplex& a)
	{
		return nvComplex(r+a.r, i+a.i);
	}
	__device__ nvComplex operator-(const nvComplex& a)
	{
		return nvComplex(r-a.r, i-a.i);
	}
	__device__ nvComplex operator-(const float real)
	{
		return nvComplex(r - real, i);
	}
	__device__ nvComplex operator/(nvComplex& a)
	{
		float len = a.squareLength();
		if (len < 1.2E-7)
			return *this;
		return nvComplex((r*a.r + i*a.i)/len, (i*a.r - r*a.i)/len);
	}
	__device__ nvComplex power3()
	{
		return nvComplex(r*r*r - 3*r*i*i, 3*r*r*i - i*i*i);
	}
	__device__ nvComplex power2()
	{
		return nvComplex(r*r - i*i, 2*r*i);
	}
};
/*==========================================================================*/

/*==========================================================================*/

__device__ float nova(int x, int y)
{
	nvComplex R(0.8f, 0.2f);
	nvComplex c(0.1f, 0);
	float jx = SCALE * (float)(DIM/2 - x)/(DIM/2);
	float jy = SCALE * (float)(DIM/2 - y)/(DIM/2);
	nvComplex z(jx, jy);
	int i;
	float delta = 0.06f;
	nvComplex oldZ(0, 0);

	for (i = 0; i < 60 && abs(delta) > 0.05f; ++i) 
	{
		z = z - R*((z.power3() - 1)/(z.power2()*3)) + c;
		delta = (z - oldZ).squareLength();
		oldZ = z;
	}

	return (1 - i/60.0f);

}

__global__ void kernel(unsigned char *ptr)
{
	// tính tọa độ từ chỉ số của block
	int y = blockIdx.y;
	int x = blockIdx.x;
	int offset = x + y * gridDim.x;
	/*
	ptr[offset*4 + 0] = 0;
	ptr[offset*4 + 1] = 128;
	ptr[offset*4 + 2] = 255*nova(x, y);
	ptr[offset*4 + 3] = 255;
	*/
	ptr[offset*4 + 0] = 
		ptr[offset*4 + 1] = 
		ptr[offset*4 + 2] = 255 * nova(x, y);
	ptr[offset*4 + 3] = 255;
}

int main( void )
{
	GlBitmap bitmap(DIM, DIM);
	unsigned char *dev_data;
	dim3 grid(DIM, DIM);

	SAFE_CALL( cudaMalloc(&dev_data, bitmap.get_size()));
	
	kernel<<<grid, 1>>>(dev_data);

	SAFE_CALL( cudaMemcpy(bitmap.get_data(), dev_data, bitmap.get_size(), cudaMemcpyDeviceToHost));
	cudaFree(dev_data);

	bitmap.display_and_exit();
	return 0;
}
#endif